Compare datasets
Camille
2024-02-15
library(raster)
library(plyr)
library(dplyr)
library(reshape)
library(ggplot2)
library(ggpmisc)
library(ggpubr)
library(viridisLite)
library(colorspace)
library(RColorBrewer)
library(sf)
library(stringr)
memory.limit(size = 160000)
## [1] Inf
#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")
col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "NA" = "#D5D8DC")1 Data
Liens utiles:
# era5vA <- read.table("era5_AI.txt") %>%
# mutate(tas = tas-273.15, pr = mpr * 60 * 60 * 24 * 365, rsds = ssrd / (60*60*24) * 1e-6 * 60*60*24) %>% # conversion from J/m2 to J/sm2 to MJ/m2/day
# select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","rsds","ea","ET0","AI","cat.AI"))
#
#
# write.table(era5vA, "era5vA.txt")
era5vA <- read.table("ERA5/era5_AI.txt") %>%
mutate(z = z, tas = t2m - 273.15, sfcWind = wind, model = "historical", source = "ERA5") %>% #pr, ET0 already in mm/year
select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn","ea","z","ET0","AI","cat.AI")) %>%
filter(lm == 1)
write.table(era5vA, "era5vA.txt")
wdcvA <- read.table("wdc_AI.txt") %>% mutate(z = z, ea = vapr, ET0 = ET0) %>% # srad in MJ/m2/day, pr and ETO already in mm/year
select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI","cat.AI")) %>%
filter(lm == 1)
write.table(wdcvA, "wdcvA.txt")
cmip6 <- read.table("cmip6_AI.txt")
cmip6vA <- cmip6 %>% subset(period == "1970_2000") %>%
group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, z) %>%
dplyr::summarise(tas = mean(tas, na.rm = T) - 273.15, pr = mean(spr, na.rm = T), # pr in mm/year
sfcWind = mean(sfcWind, na.rm = T), Rn = mean(Rn, na.rm = T),
ea = mean(ea, na.rm = T), ET0 = mean(sET0, na.rm = T), AI = mean(AI, na.rm = T)) %>% # ET0 in mm/year
ungroup() %>%
mutate(source = "CMIP6")%>%
select(c("lon","lat","Continent","Type","Name","Acronym","lm","model","source","pr","tas","sfcWind","Rn", "ea","z","ET0","AI"))
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
cmip6vA$cat.AI <- cmip6vA$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6vA$cat.AI[which(cmip6vA$ET0 < 400)] <- "Cold"
write.table(cmip6vA, "cmip6vA.txt")
era5vA <- read.table("era5vA.txt")
wdcvA <- read.table("wdcvA.txt")
cmip6vA <- read.table("cmip6vA.txt")2 Regions and Continents
colors_continents <- c("#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
"#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf",
"#aec7e8", "#ffbb78", "#98df8a", "#ff9896")
ggplot(era5vA)+geom_raster(aes(x= lon, y = lat, fill = Continent))+
scale_fill_discrete(type = colors_continents)+
theme_void()
“Continents” to remove are: * Arctic (too little gridcells) * Atlantic (ocean) * Indian (ocean) * Pacific (ocean) * Southern (ocean)
land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
plot(land_mask)
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land
elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif")
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))
ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326"))
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))
ipcc_sf <- st_as_sf(ipcc_regions)
ggplot(data = filter(ipcc_sf, Type %in% c("Land","Land-Ocean")))+
geom_sf(aes(fill=Continent))+
scale_fill_manual(values = colorvec)+
geom_sf_text(aes(label = Name), size = 3, position = position_jitter(height = 3, width = 3, seed = 1))+
#geom_sf_text(aes(label = Acronym), size = 2)+
theme_void()+
theme(legend.position = "bottom")
#ggsave("ipcc_regions.png", width = 12, height = 8, units = "cm", dpi = "retina")3 Aridity Index
comp <- merge(select(era5vA, c("lon","lat","source","Continent","z","Rn","ET0","AI","cat.AI")),
select(wdcvA, c("lon","lat","source","Continent","z","Rn","ET0","AI","cat.AI")), by = c("lon","lat")) %>%
merge(select(cmip6vA, c("lon","lat","source","Continent","z","ET0","AI","cat.AI")), by = c("lon","lat"))3.1 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000
calib <- rbind(cmip6vA, era5vA, wdcvA)
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"
ggarrange(
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
labs(title = "Wordclim, 1970-2000", fill = "")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
labs(title = "ERA5, 1970-2000", fill = "")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
labs(title = "CMIP6, 1970-2000", fill = "")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
geom_raster(aes(x = lon, y = lat), fill = "white")+
scale_fill_manual(values = col.cat)+
labs(title = "", fill = "")+
theme_void()+
theme(legend.position = "left")
)
3.2 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with country borders
calib <- rbind(cmip6vA, era5vA, wdcvA)
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"
ggarrange(
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
borders(colour= "grey20")+
labs(title = "Wordclim, 1970-2000", fill = "")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
borders(colour= "grey20")+
labs(title = "ERA5, 1970-2000", fill = "")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
borders(colour= "grey20")+
labs(title = "CMIP6, 1970-2000", fill = "")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
geom_raster(aes(x = lon, y = lat), fill = "white")+
scale_fill_manual(values = col.cat)+
labs(title = "", fill = "")+
theme_void()+
theme(legend.position = "top"),
nrows = 4, ncol = 1
)
bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <- brewer.pal(n = 11, name ="RdYlBu")
col.pr <- c("(-400,0]" = colvecrdbu[5], "(0,400]"= colvecrdbu[6],"(400,800]"= colvecrdbu[7], "(800,1.2e+03]" = colvecrdbu[8], "(1.2e+03,1.6e+03]"= colvecrdbu[9], "(1.6e+03,5e+03]" = colvecrdbu[10])
calib$pr.cat <- calib$pr %>% cut(breaks = bpr)
g.pr.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat, fill = pr.cat))+
scale_fill_manual(values= col.pr, na.translate = F)+
borders(colour= "grey20")+
labs(title = "Annual precipitation in Wordclim", fill = "mm/day")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none")
g.pr.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat, fill = pr.cat))+
scale_fill_manual(values= col.pr, na.translate = F)+
borders(colour= "grey20")+
labs(title = "Annual precipitation in ERA5", fill = "mm/day")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none")
g.pr.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat, fill = pr.cat))+
scale_fill_manual(values= col.pr, na.translate = F)+
borders(colour= "grey20")+
labs(title = "Annual precipitation CMIP6", fill = "mm/day")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none")
l.pr <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = pr.cat))+
scale_fill_manual(values= col.pr, na.translate = F)+
geom_raster(aes(x = lon, y = lat), fill = "white")+
labs(title = "", fill = "")+
theme_void()+
theme(legend.position = "top", legend.key.width = unit(2, "cm"))
#)bpr <- c(-400, 0, 400, 800, 1200, 1600, 5000)
colvecrdbu <- brewer.pal(n = 11, name ="RdYlBu")
col.et0 <- c("(-400,0]" = colvecrdbu[8], "(0,400]"= colvecrdbu[7],"(400,800]"= colvecrdbu[6], "(800,1.2e+03]" = colvecrdbu[5], "(1.2e+03,1.6e+03]"= colvecrdbu[4], "(1.6e+03,5e+03]" = colvecrdbu[3])
calib$et0.cat <- calib$ET0 %>% cut(breaks = bpr)
g.et0.wdc <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat, fill = et0.cat))+
scale_fill_manual(values= col.et0, na.translate = F)+
borders(colour= "grey20")+
labs(title = "Annual evapotranspiration Wordclim", fill = "mm/day")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none")
g.et0.era5 <- ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat, fill = et0.cat))+
scale_fill_manual(values= col.et0, na.translate = F)+
borders(colour= "grey20")+
labs(title = "Annual evapotranspiration ERA5", fill = "mm/day")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none")
g.et0.cmip6 <- ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat, fill = et0.cat))+
scale_fill_manual(values= col.et0, na.translate = F)+
borders(colour= "grey20")+
labs(title = "Annual evapotranspiration CMIP6", fill = "mm/day")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "none")
l.et0 <- ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = et0.cat))+
scale_fill_manual(values= col.et0, na.translate = F)+
geom_raster(aes(x = lon, y = lat), fill = "white")+
labs(title = "", fill = "")+
theme_void()+
theme(legend.position = "top", legend.key.width = unit(2, "cm"))ggarrange(
ggarrange(plotlist = list(g.pr.wdc,g.pr.era5, g.pr.cmip6, l.pr), nrow = 4, ncol = 1),
ggarrange(plotlist = list(g.et0.wdc,g.et0.era5, g.et0.cmip6, l.et0), nrow = 4, ncol = 1),
ncol = 2)
3.3 Comparison of WDC, ERA5 and CMIP6 for the reference period 1970-2000 with Antarctica
calib <- rbind(cmip6vA, era5vA, wdcvA)
breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")
calib$cat.AI <- calib$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
calib$cat.AI[which(calib$ET0 < 400)] <- "Cold"
ggarrange(
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
labs(title = "Wordclim, 1970-2000", fill = "")+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "ERA5")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
labs(title = "ERA5, 1970-2000", fill = "")+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "CMIP6")) + geom_raster(aes(x=lon, y = lat, fill = cat.AI))+
scale_fill_manual(values = col.cat)+
labs(title = "CMIP6, 1970-2000", fill = "")+
theme_void()+
theme(legend.position = "none"),
ggplot(data = subset(calib, source == "Worldclim")) + geom_raster(aes(x = lon, y = lat, fill = cat.AI))+
geom_raster(aes(x = lon, y = lat), fill = "white")+
scale_fill_manual(values = col.cat)+
labs(title = "", fill = "")+
theme_void()+
theme(legend.position = "left")
)
calib$same.cat <- NA
calib$cat.ref <- NA
for(i in 1:nrow(calib)){
loni = calib$lon[i]
lati = calib$lat[i]
catrefi <- subset(calib, lon == loni & lat == lati & source == "CMIP6")$cat.AI %>% as.character()
cati <- calib$cat.AI[i] %>% as.character()
calib$same.cat[i] <- ifelse(cati == catrefi, "y",
ifelse(cati!= catrefi, "n",
NA))
calib$cat.ref[i] <- catrefi
}
write.table(calib, "calib.txt")# table_calib_c <- calib %>% select(c("Continent","source","cat.AI")) %>%
# reshape2::melt(id.vars = c("Continent","source")) %>%
# group_by(Continent,source,value) %>%
# summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
# ungroup()
#
# table_calib_c$value[is.na(table_calib_c$value) == T] <- "NA"
#
# ggplot(table_calib_c)+
# geom_col(aes(x = Continent, y = perc, fill = value, position = "fill"))+
# scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
# theme_minimal()+
# theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))+
# facet_grid(rows = vars(source))+
# labs(x = "Continent", y = "%", fill = "Aridity\ncategory", y = "Proportion of aridity category")
#
table_calib <- calib %>% select(c("Continent","Name","Acronym","source","cat.AI")) %>%
reshape2::melt(id.vars = c("Continent","Name","Acronym","source")) %>%
filter(!is.na(value)) %>%
group_by(Continent,Name,Acronym,source,value) %>%
summarise(count = n()) %>% mutate(perc = round(count/sum(count)*100,1)) %>%
ungroup() %>% filter(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC"))
table_calib$value[is.na(table_calib$value) == T] <- "NA"
table_calib$Continent[which(table_calib$Continent == "EUROPE-AFRICA")] <- "EUROPE"
table_calib$Continent[which(table_calib$Continent == "CENTRAL-AMERICA")] <- "CENTRAL\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "NORTH-AMERICA")] <- "NORTH\nAMERICA"
table_calib$Continent[which(table_calib$Continent == "SOUTH-AMERICA")] <- "SOUTH\nAMERICA"
table_calib$Continent <- factor(table_calib$Continent, levels = c("AFRICA","EUROPE","ASIA","OCEANIA","NORTH\nAMERICA","CENTRAL\nAMERICA","SOUTH\nAMERICA"))
ggplot(table_calib)+
geom_col(aes(x = Acronym, y = perc, fill = value, position = "fill"))+
scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
theme_minimal(base_size = 12)+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
legend.position = "bottom",
strip.text.x = element_text(size = 6))+
facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
labs(x = "", y = "%", fill = "Aridity\ncategory")
ggplot(table_calib)+
geom_col(aes(x = Acronym, y = count, fill = value, position = "fill"))+
scale_fill_manual(values = col.cat, breaks = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","NA"))+
scale_y_continuous(position = "right")+
theme_minimal(base_size = 12)+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
legend.position = "bottom",
strip.text.x = element_text(size = 6))+
facet_grid(source ~ Continent, switch = "y", scales = "free_x", space = "free")+
labs(x = "", y = "Count", fill = "Aridity\ncategory")
library(ggrepel)
pie_calib <- table_calib %>% group_by(source, value) %>% summarise(gridcells = sum(count)) %>% mutate(percent = round(gridcells/sum(gridcells)*100, 1), lab.y = (rev(cumsum(rev(gridcells)))) - gridcells*0.5)
ggplot(pie_calib)+geom_bar(aes(x = "", y = gridcells, fill = value), width = 1, stat = "identity")+
coord_polar("y", start = 0)+
scale_fill_manual(values = col.cat, na.translate = T)+
facet_grid(cols = vars(source))+
geom_label_repel(aes(x = "", y = lab.y, label = percent), nudge_x = 0.5)+
labs(fill = "")+
theme_void()+theme(legend.position = "bottom")
4 Number of inhabitants by region
density <- raster("Population_density/gpw_v4_population_density_rev11_2020_1_deg.tif")
plot(density)
popcount <- raster("Population_density/gpw_v4_population_count_rev11_2020_1_deg.tif")
plot(popcount)
land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land
popcount.df <- popcount %>%
projectRaster(to = land_mask) %>%
as.data.frame(xy = T) %>% setNames(c("lon","lat","pop")) %>%
merge(land_mask.df, by = c("lon","lat")) %>% filter(lm == 1)
era5vApop <- merge(era5vA, popcount.df, by = c("lon","lat"))
tab.pop.by.cat <- era5vApop %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
group_by(cat.AI, Continent) %>%
summarise(pop = sum(pop))
write.csv(tab.pop.by.cat, "population.count.era5.csv")knitr::kable(tab.pop.by.cat)| cat.AI | Continent | pop |
|---|---|---|
| Arid | AFRICA | 9.407395e+07 |
| Arid | ASIA | 2.608739e+08 |
| Arid | CENTRAL-AMERICA | 1.102838e+07 |
| Arid | EUROPE | 3.523270e+05 |
| Arid | EUROPE-AFRICA | 3.354203e+07 |
| Arid | NORTH-AMERICA | 8.275538e+05 |
| Arid | OCEANIA | 2.483547e+05 |
| Arid | SOUTH-AMERICA | 2.714742e+06 |
| Cold | ASIA | 9.005723e+06 |
| Cold | EUROPE | 6.822776e+06 |
| Cold | NORTH-AMERICA | 8.555420e+05 |
| Cold | POLAR | NA |
| Cold | SOUTH-AMERICA | 6.222267e+04 |
| Dry subhumid | AFRICA | 1.131375e+08 |
| Dry subhumid | ASIA | 3.724373e+08 |
| Dry subhumid | CENTRAL-AMERICA | 2.847546e+07 |
| Dry subhumid | EUROPE | 5.195230e+06 |
| Dry subhumid | EUROPE-AFRICA | 5.933504e+07 |
| Dry subhumid | NORTH-AMERICA | 9.392606e+06 |
| Dry subhumid | OCEANIA | 7.027896e+05 |
| Dry subhumid | SOUTH-AMERICA | 1.297947e+07 |
| Humid | AFRICA | 6.854010e+08 |
| Humid | ASIA | 3.055659e+09 |
| Humid | CENTRAL-AMERICA | 1.462475e+08 |
| Humid | EUROPE | 4.999673e+08 |
| Humid | EUROPE-AFRICA | 1.452014e+08 |
| Humid | NORTH-AMERICA | 2.915502e+08 |
| Humid | OCEANIA | 1.529713e+07 |
| Humid | POLAR | 4.329376e+02 |
| Humid | SOUTH-AMERICA | 3.724919e+08 |
| Hyperarid | AFRICA | 5.328568e+07 |
| Hyperarid | ASIA | 1.970440e+07 |
| Hyperarid | EUROPE-AFRICA | 4.146658e+07 |
| Hyperarid | SOUTH-AMERICA | 5.088920e+05 |
| Semi-arid | AFRICA | 2.154724e+08 |
| Semi-arid | ASIA | 5.713649e+08 |
| Semi-arid | CENTRAL-AMERICA | 3.973955e+07 |
| Semi-arid | EUROPE | 1.841667e+06 |
| Semi-arid | EUROPE-AFRICA | 8.619193e+07 |
| Semi-arid | NORTH-AMERICA | 2.806716e+07 |
| Semi-arid | OCEANIA | 3.403064e+06 |
| Semi-arid | SOUTH-AMERICA | 2.872740e+07 |
| NA | AFRICA | 3.321493e+04 |
| NA | ASIA | 1.762974e+07 |
| NA | CENTRAL-AMERICA | 7.132541e+05 |
| NA | EUROPE | 1.481967e+05 |
| NA | EUROPE-AFRICA | 1.876606e+06 |
| NA | NORTH-AMERICA | 1.403549e+03 |
| NA | OCEANIA | 3.896168e+04 |
| NA | POLAR | NA |
| NA | SOUTH-AMERICA | 2.214226e+05 |
5 Comparison of datasets
vars <- c("lon","lat","Continent","Name","Acronym","tas","pr","sfcWind","Rn","ea","ET0","AI","cat.AI")
merge.vars <- c("lon","lat","Continent","Name","Acronym")
calib.wide <- select(era5vA, vars) %>%
merge(select(wdcvA, vars),
by = merge.vars) %>%
merge(select(cmip6vA, vars),
by = merge.vars) %>%
setNames(c(merge.vars,
"tas.era5","pr.era5", "sfcWind.era5","Rn.era5","ea.era5","ET0.era5","AI.era5","cat.AI.era5",
"tas.wdc","pr.wdc", "sfcWind.wdc","Rn.wdc","ea.wdc","ET0.wdc","AI.wdc","cat.AI.wdc",
"tas.cmip6","pr.cmip6", "sfcWind.cmip6","Rn.cmip6","ea.cmip6","ET0.cmip6","AI.cmip6","cat.AI.cmip6"))5.1 Scatterplots with Antarctica
5.1.1 ERA5 compared to WDC
Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.
5.1.2 ERA5 and WDC
ggarrange(
ggplot(calib.wide,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
knitr::kable(data.frame(precipitation = cor(calib.wide$pr.wdc, calib.wide$pr.era5, use = "pairwise.complete.obs"),
temperature = cor(calib.wide$tas.wdc, calib.wide$tas.era5, use = "pairwise.complete.obs"),
wind = cor(calib.wide$sfcWind.wdc, calib.wide$sfcWind.era5, use = "pairwise.complete.obs"),
Rn = cor(calib.wide$Rn.wdc, calib.wide$Rn.era5, use = "pairwise.complete.obs"),
ea = cor(calib.wide$ea.wdc, calib.wide$ea.era5, use = "pairwise.complete.obs"),
ET0 = cor(calib.wide$ET0.wdc, calib.wide$ET0.era5, use = "pairwise.complete.obs")),
digits = 2, caption = "Correlation coefficients between ERA5 and Worldclim") %>%
kableExtra::kable_styling(bootstrap_options = "bordered")| precipitation | temperature | wind | Rn | ea | ET0 |
|---|---|---|---|---|---|
| 0.89 | 0.99 | 0.68 | 0.67 | 1 | 0.98 |
ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
geom_density(aes(x=tas.era5, col = "ERA5"))+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none")
ggarrange(
ggplot(calib.wide) + geom_density(aes(x = tas.wdc, col = "WDC"))+
geom_density(aes(x=tas.era5, col = "ERA5"))+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
5.1.3 CMIP6 and WDC
ggarrange(
ggplot(calib.wide,aes(x = tas.wdc))+
geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = pr.wdc))+
geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = sfcWind.wdc))+
geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = Rn.wdc))+
geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ea.wdc))+
geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ET0.wdc))+
geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
5.1.4 CMIP6 and ERA5
ggarrange(
ggplot(calib.wide,aes(x = tas.era5))+
geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
scale_color_manual(values = colorvec[8])+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = pr.era5))+
geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = sfcWind.era5))+
geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = Rn.era5))+
geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ea.era5))+
geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.wide,aes(x = ET0.era5))+
geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2)
5.2 Scatterplots without Antarctica
calib.noant <- calib.wide %>% filter(lat > -60)5.2.1 ERA5 compared to WDC
Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.
5.2.2 ERA5 and WDC
ggarrange(
ggplot(calib.noant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
geom_density(aes(x=tas.era5, col = "ERA5"))+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none")
ggarrange(
ggplot(calib.noant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
geom_density(aes(x=tas.era5, col = "ERA5"))+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
5.2.3 CMIP6 and WDC
ggarrange(
ggplot(calib.noant,aes(x = tas.wdc))+
geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = pr.wdc))+
geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = sfcWind.wdc))+
geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = Rn.wdc))+
geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ea.wdc))+
geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ET0.wdc))+
geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
5.2.4 CMIP6 and ERA5
ggarrange(
ggplot(calib.noant,aes(x = tas.era5))+
geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
scale_color_manual(values = colorvec[8])+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = pr.era5))+
geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = sfcWind.era5))+
geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = Rn.era5))+
geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ea.era5))+
geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.noant,aes(x = ET0.era5))+
geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2)
5.3 Scatterplots only Antarctica
calib.ant <- calib.wide %>% filter(lat < -60)5.3.1 ERA5 compared to WDC
Rn is computed differently for ERA5 and WDC, hence the discrepancy for low values.
5.3.2 ERA5 and WDC
ggarrange(
ggplot(calib.ant,aes(x = tas.wdc))+geom_point(aes( y = tas.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, WDC", y = "Temperature in °C, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
geom_density(aes(x=tas.era5, col = "ERA5"))+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none")
ggarrange(
ggplot(calib.ant) + geom_density(aes(x = tas.wdc, col = "WDC"))+
geom_density(aes(x=tas.era5, col = "ERA5"))+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = pr.wdc))+geom_point(aes( y = pr.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = sfcWind.wdc))+geom_point(aes( y = sfcWind.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = Rn.wdc))+geom_point(aes( y = Rn.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ea.wdc))+geom_point(aes(y = ea.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ET0.wdc))+geom_point(aes(y = ET0.era5, col = "ERA5"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.era5, col = "ERA5")) +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, ERA5", col = "Dataset")+
scale_color_manual(values = c(colorvec[1], colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
5.3.3 CMIP6 and WDC
ggarrange(
ggplot(calib.ant,aes(x = tas.wdc))+
geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, WDC", y = "Temperature in °C, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = pr.wdc))+
geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right", label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, WDC", y = "Annual precipitation, mm, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = sfcWind.wdc))+
geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, WDC", y = "Wind speed m/s, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = Rn.wdc))+
geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, WDC, °C", y = "Rn - G, MJ/m2/day, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ea.wdc))+
geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, WDC, °C", y = "ea, kPa, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ET0.wdc))+
geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, WDC", y = "ET0, mm/day, CMIP6")+
scale_color_manual(values = c(colorvec[4]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2
)
5.3.4 CMIP6 and ERA5
ggarrange(
ggplot(calib.ant,aes(x = tas.era5))+
geom_point(aes(y=tas.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = tas.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Temperature in °C, ERA5", y = "Temperature in °C, CMIP6", col = "Dataset")+
scale_color_manual(values = colorvec[8])+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = pr.era5))+
geom_point(aes(y=pr.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = pr.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Annual precipitation, mm, ERA5", y = "Annual precipitation, mm, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = sfcWind.era5))+
geom_point(aes(y=sfcWind.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = sfcWind.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Wind speed m/s, ERA5", y = "Wind speed m/s, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = Rn.era5))+
geom_point(aes(y=Rn.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = Rn.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="Rn - G, MJ/m2/day, ERA5, °C", y = "Rn - G, W/MJ/m2/day, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ea.era5))+
geom_point(aes(y=ea.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ea.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ea, kPa, ERA5", y = "ea, kPa, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ggplot(calib.ant,aes(x = ET0.era5))+
geom_point(aes(y=ET0.cmip6, col = "CMIP6"), alpha =0.1, shape = 16, size = 0.5)+
stat_poly_eq(aes(y = ET0.cmip6, col = "CMIP6"), label.x = "right",label.y = "bottom") +
geom_abline(slope = 1, intercept = 0)+
labs(x="ET0, mm/day, ERA5", y = "ET0, mm/day, CMIP6", col = "Dataset")+
scale_color_manual(values = c(colorvec[8]))+
theme_minimal()+
theme(legend.position = "none"),
ncol = 3, nrow = 2)
5.4 Aridity category difference
5.4.1 ERA5 and WDC
calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.era5 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.era5, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
##
## Arid to Arid Arid to Humid
## 1359 1
## Arid to Hyperarid Arid to Semi-arid
## 18 155
## Cold to Arid Cold to Cold
## 19 10101
## Cold to Dry subhumid Cold to Humid
## 9 290
## Cold to Semi-arid Dry subhumid to Arid
## 20 23
## Dry subhumid to Cold Dry subhumid to Dry subhumid
## 5 205
## Dry subhumid to Humid Dry subhumid to Semi-arid
## 209 271
## Humid to Arid Humid to Cold
## 30 520
## Humid to Dry subhumid Humid to Humid
## 347 5436
## Humid to Semi-arid Hyperarid to Arid
## 210 230
## Hyperarid to Cold Hyperarid to Hyperarid
## 2 726
## NA to Cold Semi-arid to Arid
## 570 181
## Semi-arid to Cold Semi-arid to Dry subhumid
## 3 173
## Semi-arid to Humid Semi-arid to Semi-arid
## 71 1077
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer",
ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
ifelse(change.cat.wdc == "Arid to Cold", "dryer",
ifelse(change.cat.wdc == "Cold to Humid", "dryer",
ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
ifelse(change.cat.wdc == "Humid to arid", "wetter",
ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
ifelse(change.cat.wdc == "Humid to Cold","dryer",
ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
NA))))))))))))))))))))))))
ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
borders()+
labs(fill = "ERA5 compared to WDC", title = "Change of category")+
ylim(-55,90)+
theme_void(base_size = 15)+
theme(legend.position = "bottom"),
ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat, fill = AI.era5-AI.wdc))+
binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1),
palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")),
guide = guide_legend(label.theme = element_text(angle = 0))
)+
labs(title = "Changes in AI", fill = "ERA5 compared to WDC")+
theme_void(base_size = 15)+ylim(-55,90)+
theme(legend.position = "bottom"),
ncol = 2)
ggplot() +
geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < 277 & pr.era5 - pr.wdc > - 98), aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc > 277), aes(x=lon, y = lat), fill = "#44FF44")+
geom_raster(data = subset(calib.wide, pr.era5 - pr.wdc < - 98), aes(x=lon, y = lat), fill = "#FF4444")+
labs(title = "Most extreme differences in precipitation", subtitle = "In red, the 10% most extreme gricells in which Worldclim is wetter than ERA5; \nin green, the 10% most extreme gridcells in which ERA5 is wetter than Worldclim", fill = "precipitation in ERA5 - precipitations in WDC,\nmm/year")+
theme_void(base_size = 15)+
theme(legend.position = "bottom")
5.4.2 CMIP6 and WDC
calib.wide$same.cat.wdc <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.wdc, "y","n"))
calib.wide$change.cat.wdc <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.wdc, sep = " to "))
table(calib.wide$change.cat.wdc)
##
## Arid to Arid Arid to Dry subhumid
## 1099 9
## Arid to Hyperarid Arid to Semi-arid
## 83 150
## Cold to Arid Cold to Cold
## 25 10491
## Cold to Dry subhumid Cold to Humid
## 5 588
## Cold to Semi-arid Dry subhumid to Arid
## 23 36
## Dry subhumid to Dry subhumid Dry subhumid to Humid
## 147 196
## Dry subhumid to Hyperarid Dry subhumid to Semi-arid
## 1 350
## Humid to Arid Humid to Cold
## 77 138
## Humid to Dry subhumid Humid to Humid
## 427 5112
## Humid to Hyperarid Humid to Semi-arid
## 1 435
## Hyperarid to Arid Hyperarid to Cold
## 85 2
## Hyperarid to Hyperarid NA to Cold
## 652 570
## Semi-arid to Arid Semi-arid to Dry subhumid
## 520 146
## Semi-arid to Humid Semi-arid to Hyperarid
## 111 7
## Semi-arid to Semi-arid
## 775
calib.wide$dryer.wdc <- with(calib.wide, ifelse(change.cat.wdc == "Arid to Humid", "dryer",
ifelse(change.cat.wdc == "Arid to Semi-arid", "dryer",
ifelse(change.cat.wdc == " Arid to Dry subhumid", "dryer",
ifelse(change.cat.wdc == "Arid to Hyperarid", "wetter",
ifelse(change.cat.wdc == "Arid to Cold", "dryer",
ifelse(change.cat.wdc == "Cold to Humid", "dryer",
ifelse(change.cat.wdc == "Cold to Hyperarid", "wetter",
ifelse(change.cat.wdc == "Dry subhumid to Cold", "dryer",
ifelse(change.cat.wdc == "Dry subhumid to Humid", "dryer",
ifelse(change.cat.wdc == "Dry subhumid to Arid", "wetter",
ifelse(change.cat.wdc == "Dry subhumid to Semi-arid", "wetter",
ifelse(change.cat.wdc == "Humid to arid", "wetter",
ifelse(change.cat.wdc == "Humid to Semi-arid", "wetter",
ifelse(change.cat.wdc == "Humid to Cold","dryer",
ifelse(change.cat.wdc == "Humid to Dry subhumid","wetter",
ifelse(change.cat.wdc == "Semi-arid to Dry subhumid", "dryer",
ifelse(change.cat.wdc == "Semi-arid to Arid", "wetter",
ifelse(change.cat.wdc == "Semi-arid to Humid", "dryer",
ifelse(change.cat.wdc == "Semi-arid to Cold", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Arid", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Semi-arid", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Dry subhumid", "dryer",
ifelse(change.cat.wdc == "Hyperarid to Humid", "dryer",
NA))))))))))))))))))))))))
ggarrange(
ggplot(filter(calib.wide, same.cat.wdc == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.wdc))+
scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
borders()+
labs(fill = "CMIP6 compared to WDC", "Change of category")+
ylim(-55,90)+
theme_void(base_size = 15)+
theme(legend.position = "bottom"),
ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat, fill = AI.cmip6-AI.wdc))+
borders()+
binned_scale(aesthetics = "fill", breaks = c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")),
guide = guide_legend(label.theme = element_text(angle = 0)))+
labs(title = "Change in AI", fill = "CMIP6 compared to WDC")+
theme_void(base_size = 15)+ylim(-55,90)+
theme(legend.position = "bottom"),
ncol = 2)
5.5 CMIP6 and ERA5
calib.wide$same.cat.era5 <- with(calib.wide, ifelse(cat.AI.cmip6 == cat.AI.era5, "y","n"))
calib.wide$change.cat.era5 <- with(calib.wide, paste(cat.AI.cmip6, cat.AI.era5, sep = " to "))
table(calib.wide$change.cat.era5)
##
## Arid to Arid Arid to Dry subhumid
## 989 6
## Arid to Humid Arid to Hyperarid
## 2 230
## Arid to Semi-arid Cold to Cold
## 114 10377
## Cold to Dry subhumid Cold to Humid
## 10 745
## Dry subhumid to Arid Dry subhumid to Dry subhumid
## 14 169
## Dry subhumid to Humid Dry subhumid to Semi-arid
## 219 328
## Humid to Arid Humid to Cold
## 5 62
## Humid to Dry subhumid Humid to Humid
## 382 5456
## Humid to Semi-arid Hyperarid to Arid
## 285 14
## Hyperarid to Hyperarid NA to NA
## 725 570
## Semi-arid to Arid Semi-arid to Dry subhumid
## 511 146
## Semi-arid to Humid Semi-arid to Hyperarid
## 121 3
## Semi-arid to Semi-arid
## 778
calib.wide$dryer.era5 <- with(calib.wide, ifelse(change.cat.era5 == "Arid to Humid", "dryer",
ifelse(change.cat.era5 == "Arid to Semi-arid", "dryer",
ifelse(change.cat.era5 == " Arid to Dry subhumid", "dryer",
ifelse(change.cat.era5 == "Arid to Hyperarid", "wetter",
ifelse(change.cat.era5 == "Cold to Humid", "dryer",
ifelse(change.cat.era5 == "Dry subhumid to Cold", "dryer",
ifelse(change.cat.era5 == "Dry subhumid to Humid", "dryer",
ifelse(change.cat.era5 == "Dry subhumid to Arid", "wetter",
ifelse(change.cat.era5 == "Dry subhumid to Semi-arid", "wetter",
ifelse(change.cat.era5 == "Humid to arid", "wetter",
ifelse(change.cat.era5 == "Humid to Semi-arid", "wetter",
ifelse(change.cat.era5 == "Humid to Cold","dryer",
ifelse(change.cat.era5 == "Humid to Dry subhumid","wetter",
ifelse(change.cat.era5 == "Semi-arid to Dry subhumid", "dryer",
ifelse(change.cat.era5 == "Semi-arid to Arid", "wetter",
ifelse(change.cat.era5 == "Semi-arid to Humid", "dryer",
ifelse(change.cat.era5 == "Hyperarid to Arid", "dryer",
ifelse(change.cat.era5 == "Hyperarid to Semi-arid", "dryer",
NA)))))))))))))))))))
ggarrange(
ggplot(filter(calib.wide, same.cat.era5 == "n"))+geom_tile(aes(x = lon, y = lat, fill = dryer.era5))+
scale_fill_manual(values = c( "#ef8a62", "#67a9cf"), na.translate = F)+
borders()+
labs(fill = "CMIP6 compared to ERA5", "Change of aridity category")+
ylim(-55,90)+
theme_void()+
theme(legend.position = "bottom"),
ggplot() + geom_raster(data = calib.wide, aes(x=lon, y = lat, fill = AI.cmip6-AI.era5))+
borders()+
binned_scale(aesthetics = "fill", breaks =c(-1, -0.1, -0.01, 0, 0.01, 0.1, 1), palette = function(x) rev(c("#08519c", "#2166ac", "#67a9cf", "aliceblue","bisque", "#ef8a62", "#b2182b", "#67001f")))+
labs(title = "Change of AI", fill = "CMIP6 compared to ERA5")+
theme_void()+ylim(-55,90)+
theme(legend.position = "bottom"),
ncol = 2)
6 NDVI
# Download Gimms NDVI
library(gimms)
library(gganimate)
land_mask <- raster("Masks/land_sea_mask_1degree.nc4")
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land
#ndvi.ref <- downloadGimms(x= 1981,y= 2000,dsn = "NDVI")
list.ndvi <- list.files(path = "NDVI", pattern = "ndvi.*\\.nc4$", full.names = T)[1:2]
ndvi.rast <- rasterizeGimms(list.ndvi, keep = 0) %>% mean(na.rm = T) %>% projectRaster(land_mask) # rasterize Gimms does the quality control and removes NAs
ndvi.df <- ndvi.rast %>% as.data.frame(xy = T) %>% setNames(c("lon","lat","ndvi")) %>% merge(land_mask.df, by = c("lon","lat")) %>% filter(lm != 0)
write.table(ndvi.df, "NDVI/ndvi.df.txt")ndvi.df <- read.table("NDVI/ndvi.df.txt")
ggplot(subset(ndvi.df, lm == 1))+geom_raster(aes(x=lon, y = lat, fill = ndvi))+
scale_fill_continuous_sequential(palette = "Greens 3", rev = T)+
ylim(-55,90)+labs(fill = "NDVI")+
theme_void()+theme(legend.position = "bottom")
calibn <- merge(calib, ndvi.df, by = c("lon","lat"))
ggplot(subset(calibn, !is.na(cat.AI)))+geom_boxplot(aes(x=cat.AI, y = ndvi, col = cat.AI))+facet_grid(rows = vars(source))+
scale_color_manual(values = col.cat)+
labs(x = "", y = "NDVI")+
theme_minimal()+theme(legend.position = "none")
aov(ndvi~cat.AI, data = group_by(calibn, source)) %>% summary()
## Df Sum Sq Mean Sq F value Pr(>F)
## cat.AI 5 1026 205.17 8708 <2e-16 ***
## Residuals 43200 1018 0.02
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 23577 observations effacées parce que manquantes
range.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>%
summarise(range.ndvi = paste(round(quantile(ndvi, na.rm = T)[c(2,4)],2), collapse = " to ")) %>%
cast(source~cat.AI)
cor.ndvi <- calibn %>% filter(!is.na(cat.AI)) %>% group_by(source, cat.AI) %>%
summarise(pearson.r = round(cor(ndvi, AI, use = "pairwise.complete"),2)) %>%
cast(source~cat.AI)
knitr::kable(range.ndvi, caption = "NDVI range by source and aridity category")| source | Cold | Hyperarid | Arid | Semi-arid | Dry subhumid | Humid |
|---|---|---|---|---|---|---|
| CMIP6 | 0.3 to 0.47 | 0.09 to 0.11 | 0.1 to 0.19 | 0.17 to 0.36 | 0.24 to 0.49 | 0.4 to 0.69 |
| ERA5 | 0.27 to 0.47 | 0.09 to 0.11 | 0.11 to 0.21 | 0.19 to 0.37 | 0.24 to 0.48 | 0.42 to 0.68 |
| Worldclim | 0.29 to 0.47 | 0.09 to 0.11 | 0.1 to 0.19 | 0.18 to 0.34 | 0.28 to 0.44 | 0.45 to 0.7 |
knitr::kable(cor.ndvi, caption = "Correlation between NDVI and AI by source and aridity category")| source | Cold | Hyperarid | Arid | Semi-arid | Dry subhumid | Humid |
|---|---|---|---|---|---|---|
| CMIP6 | -0.08 | 0.15 | 0.36 | 0.38 | -0.01 | 0.31 |
| ERA5 | -0.16 | 0.08 | 0.47 | 0.30 | 0.08 | 0.38 |
| Worldclim | 0.05 | 0.05 | 0.49 | 0.39 | 0.14 | 0.52 |
knitr::knit_exit()6.0.1 Maps CMIP6 ERA5 // WDC
In red: when CMIP6 or ERA5 are warmer than Wordlclim
cowplot::plot_grid(
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = tas.era5-tas.wdc))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: ERA5 - WDC")+
theme_void()+
theme(legend.position = "bottom"),
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.wdc))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - WDC")+
theme_void()+
theme(legend.position = "bottom")
)
In blue: when CMIP6 or ERA5 are wetter than Worldclim
cowplot::plot_grid(
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = pr.era5-pr.wdc))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: ERA5 - WDC", guide = "legend")+
theme_void()+
theme(legend.position = "bottom"),
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.wdc))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - WDC", guide = "legend")+
theme_void()+
theme(legend.position = "bottom")
)
In blue: when wind speed is higher for CMIP6 or ERA5 compared to Worldclim
cowplot::plot_grid(
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = sfcWind.era5-sfcWind.wdc))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: ERA5 - WDC", guide = "legend")+
theme_void()+
theme(legend.position = "bottom"),
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.wdc))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - WDC", guide = "legend")+
theme_void()+
theme(legend.position = "bottom")
)
6.0.2 Maps CMIP6 WDC // ERA5
In red: when CMIP6 or WDC are warmer than ERA5
cowplot::plot_grid(
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = tas.wdc-tas.era5))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: WDC - ERA5")+
theme_void()+
theme(legend.position = "bottom"),
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = tas.cmip6-tas.era5))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, rev = T, name = "Temperature in °C: CMIP6 - ERA5")+
theme_void()+
theme(legend.position = "bottom")
)
In blue: when CMIP6 or WDC are wetter than ERA5
cowplot::plot_grid(
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = pr.wdc-pr.era5))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: WDC - ERA5", guide = "legend")+
theme_void()+
theme(legend.position = "bottom"),
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = pr.cmip6-pr.era5))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Precipitation in mm/y: CMIP6 - ERA5", guide = "legend")+
theme_void()+
theme(legend.position = "bottom")
)
In blue: when wind speed is higher for CMIP6 or WDC compared to ERA5
cowplot::plot_grid(
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = sfcWind.wdc-sfcWind.era5))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: WDC - ERA5", guide = "legend")+
theme_void()+
theme(legend.position = "bottom"),
ggplot(calib.wide)+
geom_raster(aes(x=lon, y = lat, fill = sfcWind.cmip6-sfcWind.era5))+
scale_fill_continuous_divergingx(palette = "RdYlBu", mid = 0, name = "Wind speed in m/s: CMIP6 - ERA5", guide = "legend")+
theme_void()+
theme(legend.position = "bottom")
)